#--------------------------------------------------------------------------
#
# Replication file for Asset specificity, Corporate Protection, and Trade Policy
# British Journal of Political Science
# Benjamin C.K. Egerod & Mogens K. Justesen
#
#--------------------------------------------------------------------------

# Preliminaries

# set working directory to the folder where the replication material is
#setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# load packages
library(sp); library(spdep);  library(splm); library(igraph)
library(stringr); library(reshape); library(reshape2); library(plyr); library(dplyr); 
library(plm); library(readr); library(ggplot2); library(grid); library(gridExtra)
library(erer); library(vcd); library(oddsratio); library(stargazer)
library(lme4); library(arm)
library(cowplot); library(broom)
library(parallel); library(ggnetwork)
library(doParallel); library(foreach); library(spatialreg)

# load firm data
firm_data <- readRDS("firm_data.rds")

# load matrix of spatial weights
w.mat2 <- readRDS("main_weights_matrix.rds")

# convert matrix to format used for analysis
panel.weight <- get.adjacency(w.mat2, attr='weight') #from edgelist to adjacency matrix
panel.weight <- as.matrix(panel.weight) # standard matrix format
panel.weight <- mat2listw(panel.weight) # from matrix to listw format

#----------------------------------------------------------------------------------
# FIGURE 1 and FIGURE 2
# Descriptive statistics
#
#----------------------------------------------------------------------------------

######################################
# FIGURE 1
# distribution of dumping decisions and duties across countries

# calculate dumping decisions
success <- firm_data %>%
  group_by(AD_CTY_NAME) %>%
  summarize(n_complaints = length(CASE_ID), # successful complaints
            n_success = sum(decision, na.rm = T)) # total n of complaints

#plot it (Panel A in full plot)
complaints1 <- ggplot(na.omit(success), aes(y = n_complaints, x = reorder(AD_CTY_NAME, n_complaints))) + 
  geom_bar(stat = "identity", colour = "black", fill = "white") +
  geom_bar(stat = "identity", aes(y = n_success, x = reorder(AD_CTY_NAME, n_complaints)),
           fill = "black") +
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"),
                     axis.text.x=element_text(angle = 60, hjust = 1)) +
  labs(x = NULL, y = "Anti-Dumping Cases  \n (log scale)") +
  coord_flip() +
  scale_y_log10(breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000))

# calculate mean of duties across countries
extent <- firm_data %>%
  group_by(AD_CTY_NAME) %>%
  summarise(extent = mean(wto_final, na.rm = T),
            n_complaints = length(CASE_ID),
            sd.extent = sd(wto_final, na.rm = T))

#Panel B: plot them with bars, using same ordering as panal A
complaints2 <- ggplot(extent, aes(y = extent, x = reorder(AD_CTY_NAME, n_complaints))) + 
  geom_bar(stat = "identity", colour = "black", fill = "black") +
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"),
                     axis.text.x=element_text(angle = 45, hjust = 1)) +
  labs(x = NULL, # remove country names from this plot
       y = "Average Duty Levied \n In Anti-Dumping Cases") +
  coord_flip() +
  scale_x_discrete(labels = NULL)

#combine plots

complaints <- grid.arrange(complaints1, complaints2, nrow = 1)

ggsave(plot = complaints, filename = "images/n_complaints.eps", 
       device = cairo_ps,
       width = 8.5, height = 6)

#####################################
# FIGURE 2
# distribution of asset specificity

#Baseline plot with homogeneous shading across entire distribution

p <- ggplot(firm_data, aes(x = mobility)) +
  geom_density(fill = "grey", alpha = .4) +
  theme_bw() + theme(panel.grid.minor = element_blank()) +
  labs(x = "Endowments of \nSpecific Assets (pct.)", y = "Density") +
  geom_vline(xintercept = mean(firm_data$mobility), lty = 2) +
  scale_x_continuous(limits = c(0,1),
                     breaks = seq(0, 1, 0.25),
                     labels = c("No Fixed \n Assets", "25 pct.",
                                "50 pct.", "75 pct.",
                                "Only Fixed \n Assets")) +
  annotate(x =.54, y = 1.25, geom = "text", label = "Mean",
           angle = 90, size = 3.5)

# add darker shading above and below 1st and 3rd quantiles
d <- ggplot_build(p)$data[[1]] # get plotted data

#calculate the two quantiles
rand1 <- quantile(firm_data$mobility, prob = 0.75)
rand2 <- quantile(firm_data$mobility, prob = 0.25)

# use new data frames to specify differently shaded areas
p <- p + geom_area(data = subset(d, x > rand1), aes(x=x, y=y), fill="grey", alpha = .5) +
  geom_area(data = subset(d, x < rand2), aes(x=x, y=y), fill="grey", alpha = .5) 

p

ggsave(plot = p, filename = "images/MobileDist.eps", device = cairo_ps,
       width = 7, height = 5)


#---------------------------------------------------------------------
#
# FIGURE 3: Plot Network of Firms Seeking Protection
#
#---------------------------------------------------------------------

gnet <- read_rds("DataForNetwork.rds")
cluster_names3 <- read_rds("LabelsForNetwork.rds")

net <- ggplot(gnet, aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_edges(color = "grey50", alpha = .2) +
  geom_nodes(size = .5) +
  theme_minimal() +
  labs(x = NULL, y = NULL) +
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  geom_text(data = cluster_names3, 
            aes(x = coord_x, y = coord_y, 
                label = industry), inherit.aes = F)

net

ggsave(plot = net, filename = "images/Network_of_Firms.eps", device = cairo_ps,
       width = 6, height = 5)

#---------------------------------------------------------------------------------
# FIGURE 4:
# Spatial Autoregressive Models
#
#---------------------------------------------------------------------------------

###########################################################################################################
# DUTY AMOUNT AS DEPENDENT VARIABLE

###############################
#bivariate models


slx.biv <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + SLX + AD_CTY_NAME,
                    data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
#summary(slx.biv)

set.seed(888)

slx.imp.biv <- impacts(slx.biv, listw = panel.weight, R = 500, zstats = T)
#summary(slx.imp.biv)

#bivariate direct effects
biv.sims <- slx.imp.biv$sres$direct
biv.sims <- as.data.frame(biv.sims)
biv.sims <- biv.sims[,1]
biv.sims <- data.frame(PE = median(biv.sims), t(quantile(biv.sims, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
biv.sims$var <- "Capital Mobility"

#bivariate indirect effects
biv.sims.indir <- slx.imp.biv$sres$indirect
biv.sims.indir <- as.data.frame(biv.sims.indir)
biv.sims.indir <- biv.sims.indir[,1]
biv.sims.indir <- data.frame(PE = median(biv.sims.indir), t(quantile(biv.sims.indir, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
biv.sims.indir$var <- "Capital Mobility"

########################
# controls

slx.mod <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + SLX + log(revenue) + log(total_assets) + log(taxation + 119067.3) 
                    + log(capital+.5)+ scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 +
                      AD_CTY_NAME+ factor(year),
                    data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
#summary(slx.mod)

set.seed(888)

slx.imp <- impacts(slx.mod, listw = panel.weight, R = 500, zstats = T)
#summary(slx.imp)

#control direct effects

cont.sims <- slx.imp$sres$direct
cont.sims <- as.data.frame(cont.sims)
cont.sims <- cont.sims[,1]
cont.sims <- data.frame(PE = median(cont.sims), t(quantile(cont.sims,probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
cont.sims$var <-c("Capital Mobility")


#control indirect effects
cont.sims.indir <- slx.imp$sres$indirect
cont.sims.indir <- as.data.frame(cont.sims.indir)
cont.sims.indir <- cont.sims.indir[,1]
cont.sims.indir <- data.frame(PE = median(cont.sims.indir), t(quantile(cont.sims.indir, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
cont.sims.indir$var <-c("Capital Mobility")

#################################
# INTRA-INDUSTRY COMPARISONS

prod2 <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + SLX + log(revenue)+
                    log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                    scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 + #SLX2 + #+ SLX3 + SLX4 + SLX5 + 
                    prod4, 
                  data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))

#summary(prod2)

set.seed(888)
prod.boots <- impacts(prod2, listw = panel.weight, R = 500, zstats = T)#summary(prod.boots)

summary(prod.boots)

# direct effects
prod.res<-as.data.frame(prod.boots$sres$direct)
prod.res <- prod.res[,1]
prod.res <- data.frame(PE = median(prod.res), t(quantile(prod.res, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
prod.res$var <- c("Capital Mobility")

#indirect effects

prod.res.indir <-as.data.frame(prod.boots$sres$indirect)
prod.res.indir <- prod.res.indir[,1:6]
prod.res.indir <- data.frame(PE = apply(prod.res.indir,2,median), t(apply(prod.res.indir,2,quantile, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
prod.res.indir$var <- c("Capital Mobility")

###################################
# gather all effect estimates

direct <- rbind(biv.sims, cont.sims, prod.res)
direct$model <- c("Bivariate", "Controls", "Product FEs")
#direct$model <- c(rep("A: Bivariate",2), rep("B: Multivariate", 10), rep("C: Product FEs", 10))
direct$type <- "Effect type: Direct"
rownames(direct)<- NULL

indirect <- rbind(biv.sims.indir, cont.sims.indir, prod.res.indir[1,])
indirect$model <- c("Bivariate", "Controls", "Product FEs")
indirect$type <- "Effect type: Indirect"
rownames(indirect) <- NULL

effects <- rbind(direct, indirect)
effects$ordering <- c(3, 2, 1, 3, 2, 1)
effects$rho.se <- c(round(slx.biv$rho.se, digits = 3), 
                    round(slx.mod$rho.se, digits = 3), 
                    round(prod2$rho.se, digits = 3), NA, NA, NA)
effects$rho <- c(round(slx.biv$rho, digits = 3), 
                 round(slx.mod$rho, digits = 3), 
                 round(prod2$rho, digits = 3), NA, NA, NA)
#effects$ordering <- c(1,2, 1, 2, 3,4,5,6,7,8,9,10, 1,2,3,4,5,6,7,8,9,10)
names(effects)[2:5] <- c("lwr95", "lwr90", "upr90", "upr95")

##############################
# plot the results

# Define some parameters for the plot
#

cols <- c("direct" = "black", "indirect" = "grey")

ann_text <- data.frame(PE = c(effects$upr95[1]+0.15,
                              effects$lwr95[2] - 0.15,
                              effects$upr95[3] + 0.15),#c(1, 0.2, 0.2, 1, 0.2, 0.2), 
                       ordering = c(3.05, 2.05, 1.05, 3, 2, 1), 
                       lab = effects$rho, lab.se = effects$rho.se,
                       type = factor(c("Effect type: Direct", "Effect type: Direct", "Effect type: Direct", 
                                       "Effect type: Indirect","Effect type: Indirect", "Effect type: Indirect"),
                                     levels = c("Effect type: Direct", "Effect type: Indirect")))
ann_text2 <- data.frame(PE = c(effects$upr95[1] +0.15,
                               effects$lwr95[2] - 0.15,
                               effects$upr95[3] + 0.15),#c(1, 0.2, 0.2, 1, 0.2, 0.2), 
                        ordering = c(2.95, 1.95, 0.95, 3, 2, 1), 
                        lab = effects$rho, lab.se = effects$rho.se,
                        type = factor(c("Effect type: Direct", "Effect type: Direct", "Effect type: Direct", 
                                        "Effect type: Indirect","Effect type: Indirect", "Effect type: Indirect"),
                                      levels = c("Effect type: Direct", "Effect type: Indirect")))

ann_text$cors <- paste(c("rho==", "rho==", "rho==", NA, NA, NA), c(ann_text$lab[1:3]),
                       sep = "")
ann_text$cors[4:6] <- NA
ann_text2$SEs <- paste("(", ann_text$lab.se, ")",
                       sep = "")
ann_text2$SEs[4:6] <- NA


#plot the  results

spat.plot <- ggplot(effects, aes(x = ordering, y = PE)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  # geom_linerange(aes(x =ordering, ymin = lwr90_2,
  #                    ymax = upr90_2),
  #                lwd = 1, position = position_dodge(width = 1/2)) +
  geom_pointrange(aes(x =  ordering, y = PE, ymin = lwr95,
                      ymax = upr95),
                  lwd = 0.65, position = position_dodge(width = 1/2),
                  shape = 21, fill = "WHITE")  +
  geom_pointrange(aes(x =  ordering, y = PE, ymin = lwr90,
                      ymax = upr90),
                  lwd = 0.95, position = position_dodge(width = 1/2),
                  shape = 21, fill = "WHITE") +
  facet_wrap(~ type, scales = "free_x") + 
  coord_flip() +
  theme_bw() + theme(panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"),
                     plot.title = element_text(hjust = 0.5)) +
  scale_colour_manual(name="Effect Type", values=cols) + 
  labs(x = NULL, y = "Marginal effect of Asset Specificity") +
  ggtitle(expression(italic("B: Modeling Duty Size"))) + 
  geom_text(data = ann_text, aes(label = cors), size = 3, parse = T) +
  geom_text(data = ann_text2, aes(label = SEs), size = 3) +
  scale_x_continuous(breaks = c(1, 2, 3),
                     labels = c("Product FEs", "Controls", "Bivariate")) 


spat.plot

#ggsave(plot = spat.plot, filename = "images/coefFe.eps", device = cairo_ps)




###########################################################################################################
# DUMPING DECISION AS DEPENDENT VARIABLE

###############################
#bivariate models


dec.biv <- lagsarlm(decision ~ I(fixed_assets/total_assets) + SLX + AD_CTY_NAME,
                    data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec.biv)

set.seed(888)

dec.imp.biv <- impacts(dec.biv, listw = panel.weight, R = 500, zstats = T)
summary(dec.imp.biv)

#bivariate direct effects
dec.biv.sims <- dec.imp.biv$sres$direct
dec.biv.sims <- as.data.frame(dec.biv.sims)
dec.biv.sims <- dec.biv.sims[,1]
dec.biv.sims <- data.frame(PE = median(dec.biv.sims), t(quantile(dec.biv.sims, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
dec.biv.sims$var <- "Capital Mobility"

#bivariate indirect effects
biv.sims.indir <- dec.imp.biv$sres$indirect
biv.sims.indir <- as.data.frame(biv.sims.indir)
biv.sims.indir <- biv.sims.indir[,1]
biv.sims.indir <- data.frame(PE = median(biv.sims.indir), t(quantile(biv.sims.indir, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
biv.sims.indir$var <- "Capital Mobility"

########################
# controls

dec.mod <- lagsarlm(decision ~ I(fixed_assets/total_assets) + SLX + log(revenue) + log(total_assets) + log(taxation + 119067.3) + log(capital+.5)+ #SLX2 + SLX3 + SLX4 + SLX5 + 
                      AD_CTY_NAME + factor(year),
                    data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec.mod)

set.seed(888)

dec.imp <- impacts(dec.mod, listw = panel.weight, R = 500, zstats = T)
summary(slx.imp)

#control direct effects

cont.sims <- dec.imp$sres$direct
cont.sims <- as.data.frame(cont.sims)
cont.sims <- cont.sims[,1:6]
cont.sims <- data.frame(PE = apply(cont.sims,2, median), t(apply(cont.sims,2,quantile,probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
cont.sims$var <-c("Capital Mobility")


#control indirect effects
cont.sims.indir <- dec.imp$sres$indirect
cont.sims.indir <- as.data.frame(cont.sims.indir)
cont.sims.indir <- cont.sims.indir[,1:6]
cont.sims.indir <- data.frame(PE = apply(cont.sims.indir,2,median), t(apply(cont.sims.indir,2,quantile, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))

cont.sims.indir$var <-c("Capital Mobility")

#################################
# INTRA-INDUSTRY COMPARISONS

dec.prod <- lagsarlm(decision ~ I(fixed_assets/total_assets) + scale(SLX) + log(revenue)+
                       log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                       scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 + #SLX2 + #+ SLX3 + SLX4 + SLX5 + 
                       prod4,
                     data = firm_data, panel.weight, method="eigen",  
                     zero.policy=TRUE, tol.solve=1.0e-13, control=list(fdHess=TRUE))
summary(dec.prod)

set.seed(888)

prod.boots <- impacts(dec.prod, listw = panel.weight, R = 500, zstats = T)
summary(prod.boots)

# direct effects
pred.res2 <-as.data.frame(prod.boots$sres$direct)
pred.res2 <- pred.res2[,1]
pred.res2 <- data.frame(PE = median(pred.res2), t(quantile(pred.res2, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
pred.res2$var <- c("Capital Mobility")

#indirect effects

pred.res2.indir <-as.data.frame(prod.boots$sres$indirect)
pred.res2.indir <- pred.res2.indir[,1:6]
pred.res2.indir <- data.frame(PE = apply(pred.res2.indir,2,median), t(apply(pred.res2.indir,2,quantile, probs=c(0.025,0.05,0.95,0.975), na.rm=TRUE)))
pred.res2.indir$var <- c("Capital Mobility")

#--------------------------------------------------
# plot dumping decision results


direct <- rbind(dec.biv.sims[1,], cont.sims[1,], pred.res2[1,])
direct$model <- c("Bivariate", "Controls", "Product FEs")
#direct$model <- c(rep("A: Bivariate",2), rep("B: Multivariate", 10), rep("C: Product FEs", 10))
direct$type <- "Effect type: Direct"

indirect <- rbind(biv.sims.indir[1,], cont.sims.indir[1,], prod.res.indir[1,])
indirect$model <- c("Bivariate", "Controls", "Product FEs")
indirect$type <- "Effect type: Indirect"

effects <- rbind(direct, indirect)
effects$ordering <- c(3, 2, 1, 3, 2, 1)
effects$rho.se <- c(round(dec.biv$rho.se, digits = 4), 
                    round(dec.mod$rho.se, digits = 4), 
                    round(dec.prod$rho.se, digits = 3), NA, NA, NA)
effects$rho <- c(round(dec.biv$rho, digits = 3), 
                 round(dec.mod$rho, digits = 3), 
                 round(dec.prod$rho, digits = 3), NA, NA, NA)
#effects$ordering <- c(1,2, 1, 2, 3,4,5,6,7,8,9,10, 1,2,3,4,5,6,7,8,9,10)
names(effects)[2:5] <- c("lwr95", "lwr90", "upr90", "upr95")



ann_text <- data.frame(PE =c(effects$upr95[1]+0.05,
                             effects$lwr95[2] - 0.05,
                             effects$upr[3] + 0.05), 
                       ordering = c(3.05, 2.05, 1.05, 3, 2, 1), 
                       lab = effects$rho, lab.se = effects$rho.se,
                       type = factor(c("Effect type: Direct", "Effect type: Direct", "Effect type: Direct", 
                                       "Effect type: Indirect","Effect type: Indirect", "Effect type: Indirect"),
                                     levels = c("Effect type: Direct", "Effect type: Indirect")))
ann_text2 <- data.frame(PE = c(effects$upr95[1]+0.05,
                               effects$lwr95[2] - 0.05,
                               effects$upr[3] + 0.05), 
                        ordering = c(2.95, 1.95, 0.95, 3, 2, 1), 
                        lab = effects$rho, lab.se = effects$rho.se,
                        type = factor(c("Effect type: Direct", "Effect type: Direct", "Effect type: Direct", 
                                        "Effect type: Indirect","Effect type: Indirect", "Effect type: Indirect"),
                                      levels = c("Effect type: Direct", "Effect type: Indirect")))

ann_text$cors <- paste(c("rho==", "rho==", "rho==", NA, NA, NA), c(ann_text$lab[1:3]),
                       sep = "")
ann_text$cors[4:6] <- NA
ann_text2$SEs <- paste("(", ann_text$lab.se, ")",
                       sep = "")
ann_text2$SEs[4:6] <- NA

#plot the  results

dec.plot <- ggplot(effects, aes(x = ordering, y = PE)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  # geom_linerange(aes(x =ordering, ymin = lwr90_2,
  #                    ymax = upr90_2),
  #                lwd = 1, position = position_dodge(width = 1/2)) +
  geom_pointrange(aes(x =  ordering, y = PE, ymin = lwr95,
                      ymax = upr95),
                  lwd = 0.65, position = position_dodge(width = 1/2),
                  shape = 21, fill = "WHITE")  +
  geom_pointrange(aes(x =  ordering, y = PE, ymin = lwr90,
                      ymax = upr90),
                  lwd = 0.95, position = position_dodge(width = 1/2),
                  shape = 21, fill = "WHITE") +
  facet_wrap(~ type, scales = "free_x") + 
  coord_flip() +
  theme_bw() + theme(panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"),
                     plot.title = element_text(hjust = 0.5)) +
  scale_colour_manual(name="Effect Type", values=cols) + 
  labs(x = NULL, y = "Marginal effect of Asset Specificity") +
  ggtitle(expression(italic("A: Modeling Dumping Decision"))) + 
  geom_text(data = ann_text, aes(label = cors), size = 3, parse = T) +
  geom_text(data = ann_text2, aes(label = SEs), size = 3) +
  scale_x_continuous(breaks = c(1, 2, 3),
                     labels = c("Product FEs", "Controls", "Bivariate"))  

dec.plot




p<-plot_grid(dec.plot, spat.plot, nrow= 2)

p

ggsave(plot = p, 
       filename = "images/results.eps", 
       device = cairo_ps,
       height = 10,
       width = 8)

##########################################
# FIGURE 5: Heterogeneous Effects

# function for computing marginal effects

marginsplotdf<-function(model, xterm, zterm, zseq){
  coefs<-coef(model)
  cov<-vcov(model)
  intterm<-ifelse(is.na(coefs[paste(xterm,zterm,sep=":")]),paste(zterm,xterm,sep=":"),paste(xterm,zterm,sep=":"))
  dy.dx<-coefs[xterm]+coefs[intterm]*zseq
  se.dy.dx<-sqrt(cov[xterm,xterm]+zseq^2*cov[intterm,intterm]+zseq*2*cov[xterm,intterm])
  CI <- data.frame(lwr95 = dy.dx-(1.96*se.dy.dx),
                   lwr90 = dy.dx-(1.68*se.dy.dx),
                   upr90 = dy.dx+(1.68*se.dy.dx),
                   upr95 = dy.dx+(1.96*se.dy.dx))
  
  margins<-data.frame(z=zseq,dydx=dy.dx,se=se.dy.dx, CI)
  return(margins)
}


#-----------------------------------------------
# Estimate interactions 

dec_rev <- lagsarlm(decision ~ mobility*log_rev+ SLX +
                      log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                      scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 +  
                      AD_CTY_NAME,
                    data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec_rev)

dec_comp <- lagsarlm(decision ~ mobility*log_comp + log_rev+ SLX +
                       log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                       scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 +  
                       AD_CTY_NAME,
                     data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec_comp)

#######
# Compute marginal effects

# marginals at levels of logged revenue
marg_rev <- marginsplotdf(dec_rev, xterm = "mobility", zterm = "log_rev",
                          zseq = seq(min(firm_data$log_rev),
                                     max(firm_data$log_rev)))

# marginals of mobility at levels of complaints against foreign firm
marg_comp1 <- marginsplotdf(dec_comp, xterm = "mobility", zterm = "log_comp",
                            zseq = seq(min(firm_data$log_comp),
                                       max(firm_data$log_comp)))

# marginals of complaints at levels of mobility
marg_comp2 <- marginsplotdf(dec_comp, xterm = "log_comp", zterm = "mobility",
                            zseq = seq(min(firm_data$mobility),
                                       max(firm_data$mobility), 0.1))


#-------------------------------------------
# plot marginal effects in Figure 5

#-----------------------------------------------
# Panel A: levels of revenue

A <- ggplot(marg_rev, aes(x = z, y = dydx))+
  geom_line()+
  geom_ribbon(aes(ymin = lwr90, ymax=upr90), fill = "grey",
              alpha=.5)+
  geom_ribbon(aes(ymin = lwr95, ymax=upr95), fill = "grey",
              alpha=.2)+
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Firm Revenue\n(Logged)", 
       y = "Marginal Effect of Asset\nSpecificity on Dumping Decision", 
       title = "A: Revenue Moderates Asset Specificity\n")

# add marginal distribution of logged revenue
a_hist <- axis_canvas(A, axis = "x") +
  geom_density(data = firm_data, aes(x = log_rev), 
               fill = "grey", alpha = .5) 
 
combined_A <- insert_xaxis_grob(A, a_hist, position = "bottom")
ggdraw(combined_A)

#-----------------------------------------------
# Panel B: levels of complaints

B <- ggplot(marg_comp1, aes(x = z, y = dydx))+
  geom_line()+
  geom_ribbon(aes(ymin = lwr90, ymax=upr90), fill = "grey",
              alpha=.5)+
  geom_ribbon(aes(ymin = lwr95, ymax=upr95), fill = "grey",
              alpha=.2)+
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "# Other Complaints Mentioning Foreign Firm\n(Logged)", 
       y = "Marginal Effect of Asset\nSpecificity on Dumping Decision", 
       title = "B: Predatory Behavior by Foreign Firm\nModerates Asset Specificity")

b_hist <- axis_canvas(B, axis = "x") +
  geom_density(data = firm_data, aes(x = log_comp), 
               fill = "grey", alpha = .5) 
# 
combined_B <- insert_xaxis_grob(B, b_hist, position = "bottom")
ggdraw(combined_B)

#-----------------------------------------------------------
# Panel C: Levels of mobility

C <- ggplot(marg_comp2, aes(x = z, y = dydx))+
  geom_line()+
  geom_ribbon(aes(ymin = lwr90, ymax=upr90), fill = "grey",
              alpha=.5)+
  geom_ribbon(aes(ymin = lwr95, ymax=upr95), fill = "grey",
              alpha=.2)+
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Firm Asset Specificity\n", 
       y = "Marginal Effect of Predatory Behavior\n on Dumping Decision", 
       title = "C: Asset Specificity Moderates\nPredatory Behavior by Foreign Firm")

c_hist <- axis_canvas(C, axis = "x") +
  geom_density(data = firm_data, aes(x = mobility), 
               fill = "grey", alpha = .5) 
# 
combined_C <- insert_xaxis_grob(C, c_hist, position = "bottom")
ggdraw(combined_C)

#--------------------------------------------------
# combine into single plot

full_plot <- plot_grid(combined_A,
                       combined_B, 
                       combined_C,
                       nrow = 1)
ggdraw(full_plot)

# export if needed
ggsave(full_plot, filename = "images/inter_mech.eps",
       device = cairo_ps,
       width = 15, height = 5)


##########################################################################
#
# ONLINE APPENDIX
#
##########################################################################


#---------------------------------------------------------------------
# APPENDIX B: Further Details on the Dataset
#---------------------------------------------------------------------

########################
# TABLE B1: Missingness Does Not Correlate with Observables

# read in data with missings
econ_dat <- readRDS("full_data.rds")

econ_dat$specific <- econ_dat$fixed_assets / econ_dat$total_assets

econ_dat$miss_spec <- ifelse(is.na(econ_dat$specific)==T, 1, 0)
econ_dat$miss_taxation <- ifelse(is.na(econ_dat$taxation)==T, 1, 0)
econ_dat$miss_revenue <- ifelse(is.na(econ_dat$revenue)==T, 1, 0)
econ_dat$miss_ass <- ifelse(is.na(econ_dat$total_assets)==T, 1, 0)



miss_mod <- lm(miss_spec ~  revenue + total_assets + taxation,
               data=econ_dat
)


miss_spec <- lm(miss_taxation ~  revenue + total_assets + specific,
                data=econ_dat
)

miss_rev <- lm(miss_revenue ~  taxation + total_assets + specific,
               data=econ_dat
)

miss_ass <- lm(miss_total_assets ~   taxation + revenue ,
               data=econ_dat
)

stargazer(miss_mod, miss_spec, miss_rev, miss_ass,
          covariate.labels = c("Revenue", "Total Assets",
                               "Taxation", "Asset Specificity"),
          omit.stat = c("f", "adj.rsq"),
          dep.var.labels = c("Missing Specificity",
                             "Missing Taxes",
                             "Missing Revenue",
                             "Missing Assets"))

#######################
# TABLE B3: Descriptive Statistics

desc_df <- dplyr::select(firm_data, decision, wto_final, mobility, total_assets, 
                         revenue, taxation, capital, SLX, SLX2, SLX3, SLX4, SLX5)
desc_df <- as.data.frame(desc_df)
stargazer(desc_df, 
          covariate.labels = c("Dumping Decision", "Duty Size", "Asset Specificity",
                               "Total Assets", "Revenue", "Taxation", "Capital",
                               "SL Specificity", "SL Revenue", "SL Total Assets",
                               "SL Tax", "SL Capital"),
          title = "Descriptive Statistics")


######################
# TABLE B4: Correlations among covariates

cor_df <- dplyr::select(firm_data, firm_id, mobility, log_rev, log_assets, log_tax, log_cap,
                        SLX, SLX2, SLX3, SLX4, SLX5)

cor_mat <- cor(cor_df[, -c(1)])
cor_mat[upper.tri(cor_mat)] <- NA

stargazer(cor_mat,
          title = "Correlations among covariates", 
          covariate.labels = c( " ", "Asset Specificity",
                                "Revenue", "Total Assets", "Taxation", "Capital",
                                "SL Specificity", "SL Revenue", "SL Assets", "SL Tax",
                                "SL Capital")
)


#------------------------------------------------------------
# APPENDIX C: The Process of Tie Formation and Spatial Autocorrelation
#-------------------------------------------------------------

#####################
# FIGURE C.1: The Distribution of Same-Good Producers Seeking Protection.

p <- ggplot(firm_data, aes(x = count_ties)) +
  geom_histogram(fill = "grey", bins = 50) +
  theme_classic() +
  labs(x = "Number of Ties", y = "Number of Firms")



ggsave(p, filename = "images/DistTies.eps", device = cairo_ps,
       width = 7, height = 5)

####################
# FIGURE C.2: The Correlates of Seeking Protection Simultaneously (A).

tie_mod <- lm(count_ties ~ mobility + log_rev +
                log_assets + log_tax + log_cap + AD_CTY_NAME,
              data = firm_data)

tidy_tie <- tidy(tie_mod)

p1 <- ggplot(tidy_tie[2:6,],
       aes(x = estimate, y = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = estimate - 1.96*std.error,
                    xmax = estimate + 1.96*std.error),
                height = 0) +
  theme_classic() +
  geom_vline(xintercept = 0, lty = 3) +
  scale_y_discrete(labels = c("Total Assets", "Capital",
                              "Revenue", "Taxation",
                              "Asset Specificty")) +
  labs(x = "Estimate", y = NULL, 
       title = "A: With Country Fixed Effects")

p1


tie_mod2 <- lm(count_ties ~ mobility + log_rev +
                log_assets + log_tax + log_cap + PRODUCT,
              data = firm_data)

tidy_tie2 <- tidy(tie_mod2)

p2 <- ggplot(tidy_tie2[2:6,],
             aes(x = estimate, y = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = estimate - 1.96*std.error,
                     xmax = estimate + 1.96*std.error),
                 height = 0) +
  theme_classic() +
  geom_vline(xintercept = 0, lty = 3) +
  scale_y_discrete(labels = c("Total Assets", "Capital",
                              "Revenue", "Taxation",
                              "Asset Specificty")) +
  labs(x = "Estimate", y = NULL, 
       title = "B: With Product Fixed Effects")

p2

p3 <- plot_grid(p1,p2)

ggsave(p3, filename = "images/CorTies.eps", device = cairo_ps,
       width = 9, height = 7.5)


####################
# FIGURE C.3: The Correlates of Seeking Protection Simultaneously (B).

tie_mod3 <- lm(count_ties ~ SLX + SLX2 + SLX3 + SLX4 + SLX5 +
                 AD_CTY_NAME,
               data = firm_data)

tidy_tie3 <- tidy(tie_mod3)

p4<- ggplot(tidy_tie3[2:6,],
            aes(x = estimate, y = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = estimate - 1.96*std.error,
                     xmax = estimate + 1.96*std.error),
                 height = 0) +
  theme_classic() +
  geom_vline(xintercept = 0, lty = 3) +
  scale_y_discrete(labels = c("SL Specificity", "SL Revenue",
                              "SL Assets", "SL Tax", "SL Capital")) +
  labs(x = "Estimate", y = NULL, 
       title = "A: With Country Fixed Effects")


tie_mod4 <- lm(count_ties ~ SLX + SLX2 + SLX3 + SLX4 + SLX5 +
                 PRODUCT,
              data = firm_data)

tidy_tie4 <- tidy(tie_mod4)

p5<- ggplot(tidy_tie4[2:6,],
       aes(x = estimate, y = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = estimate - 1.96*std.error,
                     xmax = estimate + 1.96*std.error),
                 height = 0) +
  theme_classic() +
  geom_vline(xintercept = 0, lty = 3) +
  scale_y_discrete(labels = c("SL Specificity", "SL Revenue",
                              "SL Assets", "SL Tax", "SL Capital")) +
  labs(x = "Estimate", y = NULL, 
       title = "B: With Product Fixed Effects")
p5

p6 <- plot_grid(p4,p5)

ggsave(p6, filename = "images/CorTies_slx.eps", device = cairo_ps,
       width = 9, height = 7.5)

###################
# FIGURE C.4: Diagnosing Spillover of Antidumping Duties.

I1 <- lm.morantest(lm(decision ~ mobility, data = firm_data),
            panel.weight, zero.policy = T)
I1

I2 <- lm.morantest(lm(log(wto_final+.5) ~ mobility, data = firm_data),
                   panel.weight, zero.policy = T)

I2

p7<-ggplot(firm_data, aes(x = decision, y = spat_dec)) +
  geom_point(colour = "grey") +
  geom_smooth(method = "lm", colour = "black", se = F) +
  theme_classic() +
  geom_rug() +
  labs(x = "Dumping Decision\n", 
       y = "Spatial Lag of Dumping Decision",
       title = "A: Dumping Decision") +
  annotate(geom = "text", x = 0.1, y = 50,
           label = paste(paste("Moran's I = ",
                         round(I1$statistic, digits = 1), sep = ""),
                         paste("P Value = ",
                         round(I1$p.value, digits = 1), sep = ""),
                         sep ="\n"))

p8<- ggplot(firm_data, 
       aes(x = wto_final +.5, y = spat_duty)) +
  geom_point(colour = "grey") +
  geom_smooth(method = "lm", colour = "black", se = F) +
  theme_classic() +
  geom_rug() +
  labs(x = "Duty Size\n(Log scale)", 
       y = "Spatial Lag of Duty Size\n(Variable Logged)",
       title = "B: Duty Size") +
  annotate(geom = "text", x = 3, y = 175,
           label = paste(paste("Moran's I = ",
                               round(I2$statistic, digits = 1), sep = ""),
                         paste("P Value = ",
                               round(I2$p.value, digits = 1), sep = ""),
                         sep ="\n")) +
  scale_x_log10() 

p9<-plot_grid(p7,p8)

ggsave(p9, filename = "images/MoranPlot.eps", device = cairo_ps,
       width = 12, height = 6.5)


#------------------------------------------------------------
# APPENDIX D: Details of the SAR Model Specification
#-------------------------------------------------------------

#########
# TABLE D1: SAR Coefficients

dec_rho <- c(round(dec.biv$rho, digits = 3),
             round(dec.mod$rho, digits = 3),
             round(dec.prod$rho, digits = 3))

dec_rho_se <- c(round(dec.biv$rho.se, digits = 4),
                round(dec.mod$rho.se, digits = 4),
                round(dec.prod$rho.se, digits = 4))
dec_rho/dec_rho_se

duty_rho <- c(round(slx.biv$rho, digits = 3),
              round(slx.mod$rho, digits = 3),
              round(prod2$rho, digits = 3))

duty_rho_se <- c(round(slx.biv$rho.se, digits = 4),
                 round(slx.mod$rho.se, digits = 4),
                 round(prod2$rho.se, digits = 4))

duty_rho/duty_rho_se

stargazer(dec.biv, dec.mod, dec.prod,
          slx.biv, slx.mod, prod2,
          omit = c("year", "AD", "prod", "slx", "SLX"),
          covariate.labels = c("Asset Specificity",
                               "Revenue (log)",
                               "Total Assets (log)",
                               "Taxes (log)",
                               "Capital (log)"),
          add.lines = list(c("Rho", c(dec_rho, duty_rho)),
                           c(" ", c(dec_rho_se, duty_rho_se)),
                           c("Country FE?", "Yes", "Yes", "Yes","Yes", "Yes", "Yes"),
                           c("Year FE?", "No", "Yes", "Yes", "No", "Yes", "Yes"),
                           c("Product FE?", "No", "No", "Yes", "No", "No", "Yes")
                           ),
          dep.var.labels = c("Dumping Decision", "ln Duty Size"),
          title = "SAR Coefficients",
          label = "tab:coef", no.space = T)


###########
# TABLE D2: Results Without Spatial Lags


lpm1 <- lm(decision ~ I(fixed_assets/total_assets) + AD_CTY_NAME + factor(year), 
           data = firm_data)
summary(lpm1)
c.lpm1 <- coef(lpm1)
lpm1.ci <- sqrt(diag(vcovHC(lpm1, type = "HC1")))

lpm2 <- lm(decision ~ I(fixed_assets/total_assets) + log(revenue) + log(total_assets) + 
             AD_CTY_NAME + factor(year), 
           data = firm_data)
summary(lpm2)
c.lpm2 <- coef(lpm2)[2]
lpm2.ci <- sqrt(diag(vcovHC(lpm2, type = "HC1")))

lpm3 <- lm(decision ~ I(fixed_assets/total_assets) + log(revenue) + log(total_assets)+ 
             log(taxation+ 119067.3) + log(capital+.5) + 
             AD_CTY_NAME+ factor(year), 
           data = firm_data)
summary(lpm3)
c.lpm3 <- coef(lpm3)
lpm3.ci <- sqrt(diag(vcovHC(lpm3, type = "HC1")))

###
# Size of AD duty as DV (OLS MODELS)

ols1 <- lm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + AD_CTY_NAME, data = firm_data)
summary(ols1)
c.ols1 <- coef(ols1)
ols.ci1 <- sqrt(diag(vcovHC(ols1, type = "HC1")))

ols2 <- lm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + log(revenue) + log(total_assets) + AD_CTY_NAME, 
           data = firm_data)
summary(ols2)
c.ols2<- coef(ols2)
ols.ci2 <- sqrt(diag(vcovHC(ols2, type = "HC1")))


ols3 <- lm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + log(revenue) + log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
             AD_CTY_NAME, 
           data = firm_data)
summary(ols3)
c.ols3 <- coef(ols3)
ols.ci3 <- sqrt(diag(vcovHC(ols3, type = "HC1")))

###
# set up the table of results

stargazer(lpm1, lpm2, lpm3, ols1, ols2, ols3, ci = T,
          se = list(lpm1.ci, lpm2.ci, lpm3.ci, ols.ci1, ols.ci2, ols.ci3),
          covariate.labels = c("Mobility", "Revenue (logged)", "Total Assets (logged)", "Taxation in USD", 
                               "Total Capital in USD"),
          omit = c("AD_CTY_NAME", "factor", "Constant"), df = F, no.space = T,
          omit.stat = c("ser","rsq", "f", "aic"),
          add.lines = list(c("Country FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")), 
          dep.var.labels = c("Dumped or not", "Duty (logged)"),
          title = "Capital Immobility and antidumping protection at the firm-level")





#------------------------------------------------------------
# APPENDIX E: Robustness Checks
#-------------------------------------------------------------

################
# TABLE E1: 

mod1 <- lm(scale(mobility) ~ I(lag_prot_dec*100), firm_data)
summary(mod1)


mod2 <- lm(scale(mobility) ~ I(lag_prot_dec*100) + AD_CTY_NAME + year, firm_data)
summary(mod2)


mod3 <- lagsarlm(decision ~ I(fixed_assets/total_assets) + I(lag_prot_dec*100) + SLX + log(revenue)+
                   log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) +
                   scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 + + AD_CTY_NAME +year ,
                 data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(mod3)



  mod4 <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets) +I(lag_prot_dec*100) + SLX + log(revenue)+
                   log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                   scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 +  AD_CTY_NAME + year,
                 data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(mod4)


# regression table
stargazer(mod1, mod2, mod3, mod4, 
          keep = c("fixed_assets/total_assets",
                   "lag_prot_dec"),
          covariate.labels = c("Asset Specificity",
                               "Previous Level of Protection"),
          dep.var.labels = c("Asset Specificity", "Dumping Decision", "Duty Size"),
          add.lines = list(c("Country FE?", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes"),
                           c("Firm Covariates?", "No", "No", "Yes", "Yes")),
          title = "Strategic Investment Under Expectation of Protection",
          label = "tab:sel")


##################
# TABLE E2: 


#------------------------
# regression table with foreign country dummies and foreign firm characterstics

# foreign country FEs
or_dec <- lagsarlm(decision ~ I(fixed_assets/total_assets)  + SLX + log(revenue)+
                     log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                     scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 + 
                     AD_CTY_NAME +year + INV_CTY_CODE,
                   data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(or_dec)

or_size <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + SLX + log(revenue)+
                      log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                      scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5  +
                      AD_CTY_NAME + year + INV_CTY_CODE,
                    data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(or_size)


# control for complaints against foreign firm
dec_foreign <- lagsarlm(decision ~ I(fixed_assets/total_assets) + log(av_comp)  + SLX + log(revenue)+
                          log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) +
                          scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5 
                        + AD_CTY_NAME + year,
                        data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec_foreign)



size_foreign <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets)+ log(av_comp) + SLX + log(revenue)+
                           log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                           scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5  + 
                           AD_CTY_NAME + year,
                         data = firm_data, panel.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(size_foreign)


stargazer(or_dec, or_size,
          dec_foreign, size_foreign, 
          keep = c("fixed_assets/total_assets", "av_comp"),
          covariate.labels = c("Asset Specificity",
                               "Other Complaints Re. Foreign Firm"),
          dep.var.labels = c("Dumping Decision", "Duty Size",
                             "Dumping Decision", "Duty Size"),
          add.lines = list(c("Country FE?", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes"),
                           c("Firm Covariates?", "Yes", "Yes", "Yes", "Yes")),
          title = "Characteristics of Foreign Firms and their Countries of Origin", 
          label = "tab:legit" )


##############
# FIGURE E1: The Role of Repeat Dumpers

p1 <- ggplot(firm_data, aes(x = av_comp)) +
  geom_histogram() +
  theme_classic() +
  scale_x_log10(breaks = scales::breaks_log(n = 6)) +
  annotate(geom = "text", x = 3, y = 500,
           label = paste("Median =", round(median(firm_data$av_comp,na.rm=T)))) +
  labs(x = "Number of Times a Foreign Firm is Accused of Dumping",
       y = "Count")

p1

ggsave(p1, filename = "images/DistOfComp.eps",
       device = cairo_ps, width = 6, height = 6)

p2<-ggplot(firm_data, aes(x = mobility, y = av_comp)) + 
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  theme_classic() +
  scale_y_log10(breaks = scales::breaks_log(n = 6)) +
  annotation_logticks(sides="l") + 
  geom_rug(sides = "b") +
  labs(x = "Asset Specificity",
       y = "Average Number of Complaints against\nForeign Firms mentioned in this complaint\n(Log scale)")
p2

p <- plot_grid(p1,p2)


ggsave(p2, filename = "images/SpecVsComp.eps",
       device = cairo_ps, width = 6, height = 6)

p


#############
# FIGURE E2: Repeat Complaints by Domestic Firms

n_repeat <- readRDS("RepeatComplainers.rds")

repeat_p<- ggplot(n_repeat, aes(x = n_repeat)) +
  geom_density(fill = "grey") +
  theme_classic() +
  labs(x = "Repeat Complaints from Domestic Firm\nagainst Foreign Firm",
       y = "Density")

repeat_p

ggsave(repeat_p, filename = "images/RepeatComplaint.eps",
       device = cairo_ps,
       width = 5.5, height = 3.8)


#############
# TABLE E3: 

'%!in%' <- function(x,y)!('%in%'(x,y))


w.mat2 <- readRDS("outlier_weights_matrix.rds")

excl.weight <- get.adjacency(w.mat2, attr='weight')
excl.weight <- as.matrix(excl.weight)
excl.weight <- mat2listw(excl.weight)

excl <- firm_data %>%
  dplyr::filter(AD_CTY_NAME %!in% c("Taiwan", "Philippines", 
                                    "European Union", "Japan", "Israel"))


dec_excl <- lagsarlm(decision ~ I(fixed_assets/total_assets)   + SLX + log(revenue)+
                       log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) +
                       scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5  + AD_CTY_NAME + year,
                     data = excl, excl.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec_excl)



size_excl <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + SLX + log(revenue)+
                        log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                        scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5  + AD_CTY_NAME + year,
                      data = excl, excl.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(size_excl)


stargazer(dec_excl, size_excl, 
          keep = c("fixed_assets/total_assets"),
          covariate.labels = c("Asset Specificity"),
          dep.var.labels = c("Dumping Decision", "Duty Size"),
          add.lines = list(c("Country FE?", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes"),
                           c("Firm Covariates?", "Yes", "Yes", "Yes", "Yes")),
          title = "Excluding Non-Typical Countries", 
          label = "tab:excl" )


##################
# TABLE E4: Chinese firms

sub_ch <- readRDS("china_data.rds")

ch_mod1 <- lm(decision ~ I(fixed_assets/total_assets) + log(revenue) +
                log(total_assets) + log(taxation+ 119067.3),
              data = sub_ch)

summary(ch_mod1)

ch_mod2 <- lm(log(wto_final +.5) ~ I(fixed_assets/total_assets)+ log(revenue) +
                log(total_assets) + log(taxation+ 119067.3) ,
              data = sub_ch)

summary(ch_mod2)

stargazer(ch_mod1, ch_mod2,
          covariate.labels = c("Asset Specificity", "Revenue (logged)",
                               "Total Assets (logged)", "Taxation (logged)"),
          omit.stat = c("ser","rsq", "adj.rsq", "f", "aic"),
          df = F,
          title = "Results for Chinese Firms", 
          label = "tab:china" )


##############
# FIGURE E3: Robustness to Excluding Firms

# load weights matrix again
w.mat2 <- readRDS("main_weights_matrix.rds")


cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

res <- foreach(i = 1:nrow(firm_data), 
               .combine = "rbind", 
               .packages = c("igraph", "spatialreg", "spdep")) %dopar% {
                 
                 loo <- firm_data[-c(i), ]
                 
                 loo_mat <- get.adjacency(w.mat2, attr='weight')
                 
                 loo_mat <- loo_mat[-c(i), -c(i)]
                 
                 loo.weight <- as.matrix(loo_mat)
                 loo.weight <- mat2listw(loo.weight)
                 
                 dec_loo <- lagsarlm(decision ~ I(fixed_assets/total_assets)   + SLX + log(revenue)+
                                       log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                                       scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5  + 
                                       AD_CTY_NAME + year,
                                     data = loo, loo.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
                 
                 size_loo <- lagsarlm(log(wto_final +.5) ~ I(fixed_assets/total_assets) + SLX + log(revenue)+
                                        log(total_assets) + log(taxation+ 119067.3) + log(capital+.5) + 
                                        scale_slx2 + scale_slx3 + scale_slx4 + scale_slx5  + 
                                        AD_CTY_NAME + year,
                                      data = loo, loo.weight, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
                 
                 loo_res <- c(coef(dec_loo)[3], dec_loo$rho,
                              coef(size_loo)[3], size_loo$rho) 
                 
                 # return the coefficients
                 loo_res
               }
stopCluster(cl) # shut down the cluster

res <- as.data.frame(res)

names(res) <- c("dec_loo", "dec_rho", "duty_loo", "duty_rho")

loo_res_long <- melt(res)

loo_res_long$parm <- ifelse(loo_res_long$variable == "dec_loo",
                            "A: Dumping Decision", 
                            ifelse(loo_res_long$variable == "dec_rho", "B: Rho (Decision)",
                                   ifelse(loo_res_long$variable == "duty_loo", "C: Duty Size", "D: Rho (Size)"))
)


p<-ggplot(loo_res_long, aes(x = value)) +
  geom_histogram(bins = 100) +
  facet_wrap(~parm, scales = "free_x") +
  theme_bw() + 
  labs(x = "Leave-One-Firm-Out Estimate",
       y = "Count")

p

ggsave(p, filename = "images/Loo_est.eps",
       device = cairo_ps,
       width = 5.5, height = 5.5)

#-----------------------------------------------------------------------
# APPENDIX F: Investigating the Mechanism
#-------------------------------------------------------------

w.mat2 <- readRDS("weights_alternative_mech.rds")

panel.weight2 <- get.adjacency(w.mat2, attr='weight')
panel.weight2 <- as.matrix(panel.weight2)
panel.weight2 <- mat2listw(panel.weight2)

crnt.df <- dplyr::select(econ_dat, decision, wto_final, fixed_assets, total_assets, revenue, total_assets, employees, taxation, capital,
                         year, CASE_ID, PRODUCT, AD_CTY_NAME)
firm_data2 <- as.data.frame(na.omit(crnt.df))
firm_data2 <- subset(firm_data2,log(taxation+119067.3) > 0)

firm_data2$mobility <- firm_data2$fixed_assets/firm_data2$total_assets
firm_data2$log_tax <- log(firm_data2$taxation+119067.3)
firm_data2$log_emp <- log(firm_data2$employees+.5)


############
# interactions with number of employees

dec.emp <- lagsarlm(decision ~ mobility*log_emp+ log(revenue) + 
                      log(total_assets) + log(taxation+119067.3) + 
                      log(capital +.5) + AD_CTY_NAME + factor(year),
                    data = firm_data2, panel.weight2, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec.emp)

coef_emp1 <- coef(dec.emp)[length(coef(dec.emp))]
se_emp1 <- sqrt(diag(vcov(dec.emp)))[length(coef(dec.emp))]

duty.emp <- lagsarlm(log(wto_final+.5) ~ mobility*log_emp+ log(revenue) + log(total_assets) +
                       log(taxation+119067.3) + log(capital +.5) + 
                       AD_CTY_NAME + factor(year),
                     data = firm_data2, panel.weight2, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(duty.emp)

coef_emp2 <- coef(duty.emp)[length(coef(duty.emp))]
se_emp2 <- sqrt(diag(vcov(duty.emp)))[length(coef(duty.emp))]

############
# Interactions with tax payments
dec.tax <- lagsarlm(decision ~ mobility*log_tax + log(revenue) +
                      log(total_assets) + log(capital +.5) + 
                      AD_CTY_NAME + factor(year),
                    data = firm_data2, panel.weight2, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(dec.tax)

coef_tax1 <- coef(dec.tax)[length(coef(dec.tax))]
se_tax1 <- sqrt(diag(vcov(dec.tax)))[length(coef(dec.tax))]

duty.tax <- lagsarlm(log(wto_final+.5) ~ mobility*log_tax+ log(revenue) + 
                       log(total_assets) + log(capital +.5) + 
                       AD_CTY_NAME + factor(year),
                     data = firm_data2, panel.weight2, method="eigen",  zero.policy=TRUE, tol.solve=1.0e-11, control=list(fdHess=TRUE))
summary(duty.tax)

coef_tax2 <- coef(duty.tax)[length(coef(duty.tax))]
se_tax2 <- sqrt(diag(vcov(duty.tax)))[length(coef(duty.tax))]

marg_emp_dec <- marginsplotdf(dec.emp, xterm = "mobility", zterm = "log_emp",
                          zseq = seq(min(firm_data2$log_emp),
                                     max(firm_data2$log_emp)))
marg_emp_duty <- marginsplotdf(duty.emp, xterm = "mobility", zterm = "log_emp",
                              zseq = seq(min(firm_data2$log_emp),
                                         max(firm_data2$log_emp)))
marg_tax_dec <- marginsplotdf(dec.tax, xterm = "mobility", zterm = "log_tax",
                               zseq = seq(min(firm_data2$log_tax),
                                          max(firm_data2$log_tax), 0.1))
marg_tax_duty <- marginsplotdf(duty.tax, xterm = "mobility", zterm = "log_tax",
                               zseq = seq(min(firm_data2$log_tax),
                                          max(firm_data2$log_tax), 0.1))

###########
# plot results

# taxation
tax.dec.plot <- ggplot(marg_tax_dec, aes(x = z, y = dydx)) +
  geom_line() +
  geom_ribbon(aes(ymin = (dydx - se*1.68),
                  ymax = (dydx + se*1.68)),
              alpha = .1, color = "grey") +
  geom_ribbon(aes(ymin = (dydx - se*1.96),
                  ymax = (dydx + se*1.96)),
              alpha = .2, color = "grey") +
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  geom_hline(yintercept = 0, lty = 2, col = "red") +
  geom_rug(data = firm_data2, aes(x = log_tax), inherit.aes =  F)+
  labs(y = "Effect of Asset Specificity", x = "Taxation, logged")  +
  ggtitle(expression(atop("Panel A: Effect on Dumping Decision", 
                          atop(italic("Moderated by Taxation"), "")))) +
  annotate(geom="text", x = 11.5, y = 0.5,
           label = paste(paste("Coef = ", round(coef_tax1,digits =2)),
                         paste(" SE = ", round(se_tax1,digits =2)),
                         sep ="\n"))

tax.dec.plot 

tax.duty.plot <- ggplot(marg_tax_duty, aes(x = z, y = dydx)) +
  geom_line() +
  geom_ribbon(aes(ymin = (dydx - se*1.68),
                  ymax = (dydx + se*1.68)),
              alpha = .1, color = "grey") +
  geom_ribbon(aes(ymin = (dydx - se*1.96),
                  ymax = (dydx + se*1.96)),
              alpha = .2, color = "grey") +
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  geom_hline(yintercept = 0, lty = 2, col = "red") +
  geom_rug(data = firm_data2, aes(x = log_tax), inherit.aes =  F)+
  labs(y = "Effect of Asset Specificity", x = "Taxation, logged")+
  ggtitle(expression(atop("Panel B: Effect on Duty Size", 
                          atop(italic("Moderated by Taxation"), "")))) +
  annotate(geom="text", x = 11.5, y = 2,
           label = paste(paste("Coef = ", round(coef_tax2,digits =2)),
                         paste(" SE = ", round(se_tax2,digits =2)),
                         sep ="\n"))
tax.duty.plot


# employees
emp.dec.plot <- ggplot(marg_emp_dec, aes(x = z, y = dydx)) +
  geom_line() +
  geom_ribbon(aes(ymin = (dydx - se*1.68),
                  ymax = (dydx + se*1.68)),
              alpha = .1, color = "grey") +
  geom_ribbon(aes(ymin = (dydx - se*1.96),
                  ymax = (dydx + se*1.96)),
              alpha = .2, color = "grey") +
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  geom_hline(yintercept = 0, lty = 2, col = "red") +
  geom_rug(data = firm_data2, aes(x = log_emp), inherit.aes =  F)+
  labs(y = "Effect of Asset Specificity", x = "Number of Employees, logged")+
  ggtitle(expression(atop("Panel C: Effect on Dumping Decision", 
                          atop(italic("Moderated by # of Employees"), "")))) +
  annotate(geom="text", x = 5, y = -0.1,
           label = paste(paste("Coef = ", round(coef_emp1,digits =2)),
                         paste(" SE = ", round(se_emp1,digits =2)),
                         sep ="\n"))

emp.dec.plot 

emp.duty.plot <- ggplot(marg_emp_duty, aes(x = z, y = dydx)) +
  geom_line() +
  geom_ribbon(aes(ymin = (dydx - se*1.68),
                  ymax = (dydx + se*1.68)),
              alpha = .1, color = "grey") +
  geom_ribbon(aes(ymin = (dydx - se*1.96),
                  ymax = (dydx + se*1.96)),
              alpha = .2, color = "grey") +
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  geom_hline(yintercept = 0, lty = 2, col = "red") +
  geom_rug(data = firm_data2, aes(x = log_emp), inherit.aes =  F)+
  labs(y = "Effect of Asset Specificity", x = "Number of Employees, logged")+
  ggtitle(expression(atop("Panel D: Effect on Duty Size", 
                          atop(italic("Moderated by # of Employees"), "")))) +
  annotate(geom="text", x = 5, y = -1,
           label = paste(paste("Coef = ", round(coef_emp2,digits =2)),
                         paste(" SE = ", round(se_emp2,digits =2)),
                         sep ="\n"))
emp.duty.plot


grid.arrange(tax.dec.plot, tax.duty.plot,
             emp.dec.plot, emp.duty.plot)

alt.exp <- arrangeGrob(tax.dec.plot, tax.duty.plot,
                       emp.dec.plot, emp.duty.plot)

ggdraw(alt.exp)

ggsave(plot = alt.exp, filename = "images/AlternativeExplanations.eps", device = cairo_ps,
       height = 10, width = 10)


######################
# FIGURE F2: Random Slopes by Country


cty.slope <- lmer(log(wto_final+1) ~ mobility + log_rev + log_assets + log_tax + log_cap + AD_CTY_NAME + (1 + mobility| AD_CTY_NAME) , data = firm_data)
#summary(cty.slope)

fx <- coef(cty.slope)$AD_CTY_NAME

cty.fx <- data.frame(cty = rownames(fx), estimate= fx[,2], spec = "DV: Duty Size")

bayes.mod <- sim(cty.slope)
bayes.fx <- bayes.mod@ranef$AD_CTY_NAME
bayes.fx <- bayes.fx[,,2]

CI.func <- function(x){
  CI <- quantile(x, probs = c(0.025, 0.975))
  return(CI)
}

cty.fx <- data.frame(cty = colnames(bayes.fx), t(apply(bayes.fx, 2, FUN = CI.func)),
                     estimate = fx[,2], spec = "DV: Duty Size")


# binary DV
cty.slope.bin <- lmer(decision ~ mobility + log_rev + log_assets + log_tax + log_cap + AD_CTY_NAME + (1 + mobility| AD_CTY_NAME) , data = firm_data)
#summary(cty.slope.bin)
fx <- coef(cty.slope.bin)$AD_CTY_NAME

binary.fx <- data.frame(cty = rownames(fx), estimate = fx[,2], spec = "DV: Decision")

bayes.mod <- sim(cty.slope.bin)
bayes.fx <- bayes.mod@ranef$AD_CTY_NAME
bayes.fx <- bayes.fx[,,2]

bin.fx <- data.frame(cty = colnames(bayes.fx), t(apply(bayes.fx, 2, FUN = CI.func)),
                     estimate = binary.fx[,2], spec = "DV: Decision")

#combining and plotting
cty.fx <- rbind(cty.fx, bin.fx)

segm.dec <- data.frame(x=0.31,
                       spec = factor(c("DV: Decision", "DV: Duty Size"),
                                     levels = c("DV: Decision", "DV: Duty Size")))

cty.plot <- ggplot(cty.fx, aes(x = estimate, y = reorder(cty, estimate))) +
  geom_point() +
  geom_line() +
  theme_bw() +
  geom_errorbarh(aes(xmin = X2.5., xmax = X97.5.), height = 0) +
  facet_wrap(~ spec) +
  labs(x= "Country-specific coefficient on Asset Specificity", y = NULL) +
  geom_vline(xintercept = 0, lty = 2)

cty.plot

ggsave(plot = cty.plot, 
       filename = "images/country_slopes.eps",
       device = cairo_ps,
       width = 6,
       height = 6)


